Carto-leaflet

Author

Claude Grasland

library(knitr)
## Global options
options(max.print="80")
opts_chunk$set(echo=TRUE,
               cache=FALSE,
               prompt=FALSE,
               tidy=FALSE,
               comment=NA,
               message=FALSE,
               warning=FALSE,
               options(scipen=999))
opts_knit$set(width=75)

# Packages utilitaires
library(dplyr)

Attaching package: 'dplyr'
The following objects are masked from 'package:stats':

    filter, lag
The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union
library(rmdformats)

# Packages graphiques
library(ggplot2)
library(RColorBrewer)

#packages cartographiques 
library(sf)
Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
library(leaflet)
library(htmlwidgets)
library(htmltools)

# Appel d'API
library(httr)
library(jsonlite)
library(geojsonsf)

Premiers pas

OBJECTIFS : Ce cours propose de fournir les bases élémentaires du logiciel Leaflet. Il est très largement inspiré d’un article d’Elena Salette publié sur l’excellent site de formation ThinkR et intitulé Cartographie interactive : comment visualiser mes données spatiales de manière dynamique avec leaflet ?

BUG WARNING : Il peut arriver que la transformation du fichier .Rmd en .html ne s’opère pas et que vous voyiez apparaître le message d’erreur suivant RMarkdown cannot knit: html_dependency not found. Ce message d’erreur persiste même après avoir corrigé le code… ce qui est très pénible. Après avoir cherché sur les forums de discussion, j’ai trouvé une première réponse sur stackoverflow qui consiste simplement à aller sur la flèche descendnate à droite du bouton knitr et effectuer un clear knitr cache avant de relancer le Knitr. Apparemment ça marche, sans que je sache bien pourquoi. Mais la solution la plus efficace semble être d’insérer une option cache=FALSE dans les options globales du fichier Markdown. Cela va sans doute un peu ralentir l’affichage des pages HTML,mais évite les problèmes. On pourra toujours rétablir cache=TRUE si nécessaire

Notre premier objectif très limité sera de construire une carte interactive utilisant le fonds de carte OpenStreetMap que l’on pourra zoomer en avant ou en arrière. La carte comportera la localisation de la place de la gare à Sucy-en-Brie avec une “épingle” de localisation comportant une photographie de la gare et un petit texte de promotion de celle-ci.

Lancement avec leaflet()

Nous allons avoir besoin des packages suivants :

  • leaflet puisque c’est l’objet même du cours !
  • dplyr afin de pouvoir construire des programmes utilisant des pipes %>%
  • sf pour charger des fonds de carte de différents types (points, lignes polygones)
  • htmltools et htmlwidgets pour ajouter des popups interactifs sur notre carte

Pour vérifier que le package leaflet est bien installé, nous créons une première carte (vide !)

map <- leaflet()

map

Et il n’y a … RIEN ! si ce n’est un bouton de zoom

Remplissage avec addTiles()

On ajoute sur ce fond de carte vide des “tuiles” cartographiques qui sont des images se modifiant selon l’échelle pour apporter plus ou moins de détails. Par défaut, le fonds de carte de référence est le fonds OpenStreetMap

library(leaflet)

map <- leaflet() %>%
          addTiles()

map

La carte est désormais interactive et on peut effectuer des zooms ou se déplacer.

Calage avec setView()

Nous allons ensuite choisir un point de référence, par exemple la place de la gare à Sucy-en-Brie. Pour trouver les coordonnées de latitude et longitude, la solution la plus simple est d’utiliser Google Maps puis de zoomer sur la zone d’étude et enfin d’effectuer un click droit avec la souris sur le point dont on cherche les coordonnées pour obtenir dans un popup les coordonnées recherchées :

coordonnnées de la place de la gare de Sucy On peut alors procéder à une double opération de centrage de notre carte et de définition d’une échelle d’observation afin que la carte produite par leafletcouvre bien la zone qui nous intéresse. Cette double opération est réalisée à l’aide de la fonction setView() assortie des trois paramètre suivants :

  • lng = pour la longitude
  • lat = pour la latitude
  • zoom = pour le degré d’aggrandissement de la carte de 1 pour le Monde entier à 20 pour une vision ulra locale
map <- leaflet() %>% 
          addTiles() %>%
          setView(lat = 48.77141, lng=2.50870, zoom = 17)

map

Une fois qu’on a vérifié le centrage avec un zoom fort (ici 17), on peut refaire la carte en utilisant un zoom plus faible, par exemple un zoom de 12 permettant de visualiser toute la commune de Sucy et les communes voisines.

map <- leaflet() %>% 
          addTiles() %>%
          setView(lat = 48.77141, lng=2.50870, zoom = 12)

map

Personalisation avec addProviderTiles()

Les tuiles OpenStreetMap qui servent de fonds de carte par défaut peuvent être remplacés par des tuiles personalisées fournies par des producteurs publics ou privés. On peut obtenir la liste des tuiles disponibles en tapant providers dans la console de R studio et les tester une par une. Mais il est souvent plus simple et plus rapide d’aller visualiser les tuiles disponibles sur ce site web où l’on peut centrer le monde sur sa zone d’étude et voir ce que donnent les différentes familles de tuiles.

A titre d’exemple, les tuiles OpenTopoMap permettent de voir la carte topographique de type IGN en couleur :

map <- leaflet() %>% 
            addProviderTiles('OpenTopoMap') %>%
          setView(lat = 48.77141, lng=2.50870, zoom = 12)

map

La couche Esri.WorldTopoMap fournit également une imagerie précise mais de couleurs plus neutre que les tuiles OpenTopoMap ou OpenStreetMap , ce qui sera intéressant si on superspose des marqueurs de couleur vive.

map <- leaflet() %>% 
            addProviderTiles('Esri.WorldTopoMap') %>%
          setView(lat = 48.77141, lng=2.50870, zoom = 12)
map

Affichage d’un point avec addMarkers()

L’usage le plus fréquent de leafletconsiste à ajouter des éléments de localisation ponctuelle appelés markerset de rendre ces objets ponctuels interactifs avec l’ouverture de fenêtres popupslorsqu’on clique dessus avec la souris. On va donc voir pas à pas comment construire de telles cartes interactives en partant du cas le plus simple (marqueur unique) pour aller vers les cas plus complexes (ensemble de marqueurs de taille, couleur et formes différentes).

Nous allons commencer par indiquer l’emplacement de la place de la gare de Sucy-en-Brie sur notre carte précédente à l’aide de la fonction addMarkers() :

map <- leaflet() %>% 
            addProviderTiles('Esri.WorldTopoMap') %>%
            setView(lat = 48.77141, lng=2.50870, zoom = 12) %>% 
            addMarkers(lat = 48.77141, lng=2.50870)
map

On constate que le marqueur donne bien la position choisi mais n’est pas interactif. Il faut ajouter plus de paramètres pour assurer l’interactivité.

Ajout d’un labelou d’un popup

On peut définir deux comportements d’un marker selon que la souris ne fait que passer dessus (label) ou selon que l’utilisateur effectue un click sur marker et déclenche l’ouverture d’une fenêtre (popup). Dans sa version la plus simple, l’interactivité consiste à ajouter une chaîne de caractère à ces deux paramètres.

icone_gare <-makeIcon(iconUrl = "img/gare_sucy_coord_googlemap.png")
map <- leaflet() %>% 
            addProviderTiles('Esri.WorldTopoMap') %>%
            setView(lat = 48.77141, lng=2.50870, zoom = 12) %>% 
            addMarkers(lat = 48.77141, lng=2.50870,
                      # En passant la souris
                      label = "GARE DE SUCY-BONNEUIL", 
                      # En cliquant sur l'icone
                       popup = "La gare RER A de Sucy Bonneuil est bien reliée aux communes 
                                 environnantes par un réseau de bus partant dans toutes les directions")
map

Amélioration du popup

Mais on peut faire beaucoup mieux, notamment pour la fenêtre popupqui peut prendre la forme d’une mini-page web dont on fixe le contenu en html avec la fonction paste0() et les dimensions avec le sous-paramètre popupOptions().

# Préparation de la fenêtre Popup
    my_popup = paste0(
      "<b> LA GARE DE SUCY",
      "</b><br/><img src=https://upload.wikimedia.org/wikipedia/commons/thumb/6/68/Gare_Sucy-en-Brie.jpg/1200px-Gare_Sucy-en-Brie.jpg width='200px'><br/>",
      "La gare RER A de Sucy Bonneuil est bien reliée aux communes 
                                 environnantes par un réseau de bus partant dans toutes les directions.")


  
# Réalisation de la carte
map <- leaflet() %>% 
            addProviderTiles('Esri.WorldTopoMap') %>%
            setView(lat = 48.77141, lng=2.50870, zoom = 12) %>% 
            addMarkers(lat = 48.77141, lng=2.50870,
                      # En passant la souris
                      label = "GARE DE SUCY-BONNEUIL", 
                      # En cliquant sur l'icone
                       popup = my_popup, 
                      # Quelques options de la popup
                        popupOptions = 
                      list(maxHeight = 150, maxWidth = 200))
map

Prolongements

Et voila, le tour est joué. Il faut maintenant réfléchir à la façon de construire une carte comportant un ensemble d’épingles similaires avec des couleurs ou des formes différentes, des messages différents, des photographies variées … Il ne sera alors évidemment pas possible d’ajouter une commande addMarkers() pour chaque épingle si la carte en comporte des centaines.

Si vous avez bien compris ce cours, vous pourrez trouver des réponses en lisant de façon autonome le reste de l’article dont nous nous somme inspiré : Cartographie interactive : comment visualiser mes données spatiales de manière dynamique avec leaflet ?

Cartographie de points

Nous allons prendre comme exemple la base des équipements sportifs qui sera utilisée par la suite pour le projet de fin de semestre. Cette base peut être obtenue à travers deux API différentes :

  • Base permanente des équipements sportifs et de loisirs : il s’agit d’une base qui est mise à jour en temps réel et est donc très utile lorsque l’on veut actualiser sa source régulièrement : https://equipements.sports.gouv.fr/explore/dataset/data-es/table/

  • Base des équipements sportifs et de loisirs en 2017 : il s’agit d’une photographie de la situation en 2017 qui ne fait plus l’objet de maintenance. Elle est accessible par le site public.opendasoft qui a déjà été vue dans le premier cours https://public.opendatasoft.com/explore/dataset/res_equipements_2017/information/?disjunctive.caracteristiques&sort=equdatecreation

Préparation des données

On télécharge à titre d’exemple la base des équipements sportifs et de loisirs de la commune de Sucy-en-Brie en 2017 à l’aide du site public.opendasoft. On utilise le format d’importation geojson et on crée un objet de type spatial features. On ne sélectionne que quelques caractéristiques des équipements

myurl <- "https://public.opendatasoft.com/api/explore/v2.1/catalog/datasets/res_equipements_2017/exports/geojson?lang=fr&refine=comlib%3A%22Sucy-en-Brie%22&timezone=Europe%2FBerlin"

don<-geojson_sf(myurl) %>% st_as_sf()

sel<-don %>% select(id = insnumeroinstall,
                    ins = insnom,
                    equ = equnom,
                    typ =  equipementtypelib,
                    anc = equanneeservice,
                    sup = equsurfaceevolution
                    ) %>%
              mutate(anc=as.numeric(anc),
                     sup=as.numeric(sup)) %>%
              filter(is.na(anc)==F,
                     is.na(sup)==F)


kable(head(sel))
id ins equ typ anc sup geometry
940710003 Aire de Proximite des Noyers City Stade Plateau EPS/Multisports/city-stades 1999 364 POINT (2.50918 48.77929)
940710018 Parc des Sports de Sucy en Brie 2 Courts de Tennis Exterieurs 1 et 2 Court de tennis 1969 648 POINT (2.53637 48.75876)
940710002 Stade Paul Meyer 2 Terrains de Football D et E Eclaires Terrain de football 1980 6825 POINT (2.50398 48.77557)
940710002 Stade Paul Meyer Terrain de Rugby Honneur A Non Eclaire Terrain de rugby 1980 8400 POINT (2.50548 48.77563)
940710016 Gymnase du Piple Le Gymnase Salle multisports 1996 800 POINT (2.52598 48.75997)
940710018 Parc des Sports de Sucy en Brie Terrain de Rugby 2 Terrain de rugby 1988 7080 POINT (2.54191 48.7595)

On décide de repérer uniquement les terrains de rugby, de football et de tennis

sel <- sel %>% filter(typ %in% c("Court de tennis","Terrain de football","Terrain de rugby")) %>%
              mutate(typ=as.factor(typ))
summary(sel)
      id                ins                equ           
 Length:27          Length:27          Length:27         
 Class :character   Class :character   Class :character  
 Mode  :character   Mode  :character   Mode  :character  
                                                         
                                                         
                                                         
                  typ          anc            sup                geometry 
 Court de tennis    :11   Min.   :1969   Min.   : 640   POINT        :27  
 Terrain de football:11   1st Qu.:1980   1st Qu.: 648   epsg:4326    : 0  
 Terrain de rugby   : 5   Median :1980   Median :5225   +proj=long...: 0  
                          Mean   :1983   Mean   :4166                     
                          3rd Qu.:1980   3rd Qu.:7005                     
                          Max.   :2016   Max.   :8625                     

Cartographie des localisations

On commence par créer une carte des localisations des installations sportives avec AddCircleMarkers()

# Réalisation de la carte
map <- leaflet() %>% 
            addProviderTiles('Esri.WorldTopoMap') %>%
            setView(lat = 48.77, lng=2.53, zoom = 13) %>%
             addCircleMarkers(data=sel)

map

Réglage de la taille des cercles

On règle la taille des cercles en fonction de la surface des installations

# Calcul du diamètre des cercles
  sel$myradius <-10*sqrt(sel$sup/max(sel$sup,na.rm=T))
  
# Réalisation de la carte
map <- leaflet() %>% 
            addProviderTiles('Esri.WorldTopoMap') %>%
            setView(lat = 48.77, lng=2.53, zoom = 13) %>%
  
             addCircleMarkers(data=sel,
                              radius= ~myradius,    # diamètre
                              stroke=FALSE,         # pas de bordure           
                              fillOpacity = 0.5)    # opacité 
            
                              

map

Réglage de la couleur des cercles

On fait varier la couleur des cercles en fonction due l’ancienneté des installations

# Calcul du diamètre des cercles
  sel$myradius <-10*sqrt(sel$sup/max(sel$sup,na.rm=T))

# Choix des classes 
    mycut<-c(1960,1970,1980,2000,2010,2016)
    
# Choix de la palette (c'est une fonction !)
   mypal <- colorBin('Spectral', 
                       reverse = T,
                       sel$anc,
                       bins=mycut)
  
# Réalisation de la carte
map <- leaflet() %>% 
            addProviderTiles('Esri.WorldTopoMap') %>%
            setView(lat = 48.77, lng=2.53, zoom = 13) %>%
  
             addCircleMarkers(data=sel,
                              radius= ~myradius,    # diamètre
                              stroke=TRUE,          # bordure   
                              weight=1  ,           # épaisseur de la bordure
                              color= "black",      # couleur de la bordure
                              opacity = 0.7  ,       # opacité de la bordure 
                              fillOpacity = 0.5,    # opacité 
                              fillColor = ~mypal(anc)
                              )    %>%
              addLegend(data = sel,
                      pal = mypal, 
                      title = "Ancienneté",
                      values =~anc, 
                      position = 'topright') 
            
                              

map

Ajout d’un popup d’information

On rajoute un popup pour afficher toutes les informations sur chaque terrain

# Calcul du diamètre des cercles
  sel$myradius <-10*sqrt(sel$sup/max(sel$sup,na.rm=T))

# Choix des classes 
    mycut<-c(1960,1970,1980,2000,2010,2016)
    
# Choix de la palette (c'est une fonction !)
   mypal <- colorBin('Spectral', 
                       reverse = T,
                       sel$anc,
                       bins=mycut)
# Préparation des popups
      mypopups <- lapply(seq(nrow(sel)), function(i) {
      paste0(  paste("Installation: " ,sel$ins[i]), '<br>', 
               paste("Equipement  : " ,sel$equ[i]), '<br>', 
               paste("Type        : " ,sel$typ[i]), '<br>',               
               paste("Surface     : ", sel$sup[i]), '<br>',
               paste("Ancienneté  : " ,sel$anc[i]))
            
            })
      mypopups<-lapply(mypopups, htmltools::HTML)  
      
  
# Réalisation de la carte
map <- leaflet() %>% 
            addProviderTiles('Esri.WorldTopoMap') %>%
            setView(lat = 48.77, lng=2.53, zoom = 13) %>%
  
             addCircleMarkers(data=sel,
                              radius= ~myradius,       # diamètre
                              stroke=TRUE,             # bordure   
                              weight=1  ,              # épaisseur de la bordure
                              color= "black",          # couleur de la bordure
                              opacity = 0.7  ,         # opacité de la bordure 
                              fillOpacity = 0.5,       # opacité du remplissage
                              fillColor = ~mypal(anc), # couleur de remplissage
                               popup = mypopups,       # Popup !
                              )    %>%
              addLegend(data = sel,
                      pal = mypal, 
                      title = "Ancienneté",
                      values =~anc, 
                      position = 'topright') 
            
                              

map

Choix des tuiles

On fait varier les tuiles pour offrir la possibilité de visualiser la position des maisons sur un plan ou sur une photo aérienne.

# Calcul du diamètre des cercles
  sel$myradius <-10*sqrt(sel$sup/max(sel$sup,na.rm=T))

# Choix des classes 
    mycut<-c(1960,1970,1980,2000,2010,2016)
    
# Choix de la palette (c'est une fonction !)
   mypal <- colorBin('Spectral', 
                       reverse = T,
                       sel$anc,
                       bins=mycut)
# Préparation des popups
      mypopups <- lapply(seq(nrow(sel)), function(i) {
      paste0(  paste("Installation: " ,sel$ins[i]), '<br>', 
               paste("Equipement  : " ,sel$equ[i]), '<br>', 
               paste("Type        : " ,sel$typ[i]), '<br>',               
               paste("Surface     : ", sel$sup[i]), '<br>',
               paste("Ancienneté  : " ,sel$anc[i]))
            
            })
      mypopups<-lapply(mypopups, htmltools::HTML)  
      
  
# Réalisation de la carte
map <- leaflet() %>% 
               # Tuiles
               addTiles(group = "OSM ") %>%
               addProviderTiles('Esri.WorldTopoMap', group = "ESRI topo.") %>%
               addProviderTiles('Esri.WorldImagery', group = "ESRI photo.") %>%
              # Contrôle des tuiles
               addLayersControl( baseGroups = c("OSM","ESRI topo.","ESRI photo."),
                                 position = "bottomright") %>%
            setView(lat = 48.77, lng=2.53, zoom = 13) %>%
  
             addCircleMarkers(data=sel,
                              radius= ~myradius,       # diamètre
                              stroke=TRUE,             # bordure   
                              weight=1  ,              # épaisseur de la bordure
                              color= "black",          # couleur de la bordure
                              opacity = 0.7  ,         # opacité de la bordure 
                              fillOpacity = 0.5,       # opacité du remplissage
                              fillColor = ~mypal(anc), # couleur de remplissage
                               popup = mypopups,       # Popup !
                              )    %>%
              addLegend(data = sel,
                      pal = mypal, 
                      title = "Ancienneté",
                      values =~anc, 
                      position = 'topright') 
            
                              

map

Ajout d’une couche de polygones

On décide d’ajouter le polygone donnant le contour de la commune de Sucy-en-Brie ainsi que celui des quarties Iris.

com<-readRDS("data/RP/map_com.RDS") %>% filter(com_name=="Sucy-en-Brie") 

iris <- readRDS("data/RP/map_iris.RDS") %>% filter(com_name=="Sucy-en-Brie")

# Calcul du diamètre des cercles
  sel$myradius <-10*sqrt(sel$sup/max(sel$sup,na.rm=T))

# Choix des classes 
    mycut<-c(1960,1970,1980,2000,2010,2016)
    
# Choix de la palette (c'est une fonction !)
   mypal <- colorBin('Spectral', 
                       reverse = T,
                       sel$anc,
                       bins=mycut)
# Préparation des popups
      mypopups <- lapply(seq(nrow(sel)), function(i) {
      paste0(  paste("Installation: " ,sel$ins[i]), '<br>', 
               paste("Equipement  : " ,sel$equ[i]), '<br>', 
               paste("Type        : " ,sel$typ[i]), '<br>',               
               paste("Surface     : ", sel$sup[i]), '<br>',
               paste("Ancienneté  : " ,sel$anc[i]))
            
            })
      mypopups<-lapply(mypopups, htmltools::HTML)  
      
  
# Réalisation de la carte
map <- leaflet() %>% 
               # Tuiles
               addTiles(group = "OSM ") %>%
               addProviderTiles('Esri.WorldTopoMap', group = "ESRI topo.") %>%
               addProviderTiles('Esri.WorldImagery', group = "ESRI photo.") %>%
              # Contrôle des tuiles
               addLayersControl( baseGroups = c("OSM","ESRI topo.","ESRI photo."),
                                 position = "bottomright") %>%
            setView(lat = 48.77, lng=2.53, zoom = 13) %>%
              addPolygons(data = iris,
                          fillOpacity = 0.1,
                          color="red",
                          weight = 0.5) %>%
               addPolygons(data = com,
                          fillOpacity = 0,
                          color="red",
                          weight = 2) %>%
  
             addCircleMarkers(data=sel,
                              radius= ~myradius,       # diamètre
                              stroke=TRUE,             # bordure   
                              weight=1  ,              # épaisseur de la bordure
                              color= "black",          # couleur de la bordure
                              opacity = 0.7  ,         # opacité de la bordure 
                              fillOpacity = 0.5,       # opacité du remplissage
                              fillColor = ~mypal(anc), # couleur de remplissage
                               popup = mypopups,       # Popup !
                              )    %>%
              addLegend(data = sel,
                      pal = mypal, 
                      title = "Ancienneté",
                      values =~anc, 
                      position = 'topright') 
            
                              

map

Exercice d’application

Construire une carte interactive de la localisation des courts de tennis en Ile de France en utilisant les données les plus récentes.

Acquisition des données

Importez un fichier de tous les courts de tennisde France au format geojson que vous transformerez au format sf. Puis sélectionnez les quatre départements 75, 92, 93, 94

id com ins equ typ anc sup geometry
920040011 92004 Tennis Club du Menil Court de Tennis Couvert Court de tennis 1983 648 POINT (2.28549 48.92151)
920250007 92025 Tennis Club Amiot Mini Tennis Court de tennis 2004 93 POINT (2.24056 48.92773)
920400006 92040 Tennis de Issy-Les-Moulineaux Courts de Tennis N° 2, 3, 4, 5 Court de tennis 1987 608 POINT (2.25639 48.81656)
920400047 92040 Halle des Sports Christiane Guillaume Terrains de Tennis Decouverts Court de tennis 2007 594 POINT (2.26426 48.82562)
930620002 93062 Centre Boulevard du Nord Courts de Tennis 4 Court de tennis 1990 648 POINT (2.51619 48.90288)
930630003 93063 Complexe Sportif Paul Baldit Mini Tennis Court de tennis 1966 135 POINT (2.44047 48.88761)

Première carte : localisation

Construisez une première carte indiquant simplement la localisation des courts de tennis sous la forme de cercles rouge de rayon 4 pixels et d’opacité 0.5. La carte devra être centrée sur Paris et couvrir Paris et la Petite Couronne.

Deuxième carte : Ancienneté

Construisez une seconde carte permettant de repérer l’ancienneté des terrains de tennis.